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■ The short-term memory effects recently observed in vibration-induced compaction of granular 

materials are studied, ft is shown that they can be explained by means of quite plausible hypothesis 
about the mesoscopic description of the evolution of the system. The existence of a critical time 
separating regimes of "anomalous" and "normal" responses is predicted. A simple model fitting into 
the general framework is analyzed in the detail. The relationship between this work and previous 
studies is discussed. 
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I. INTRODUCTION 

Experiments have shown that when a loose packing of grains is submitted to vertical vibration or "tapping" slowly 
approaches a steady state of higher packing fraction QJ2J]. The final steady density is a decreasing function of the 
dimcnsionlcss parameter characterizing the vibration intensity. Moreover, the relaxation is slower for smaller vibration 
intensity. In the time evolution of the system neither convection effects nor oscillatory behaviour are observed. The 
study of the kinetics of compaction is important both from a formal point of view and because of its economical 
relevance in many industrial processes. Most of the peculiar behaviours exhibited by granular materials submitted 
' to vibration or tapping processes show a great similarity with conventional structural glasses. This includes slow 
relaxation, annealing properties, and hysteresis effects. 

The first study of the response of a granular system to a sudden change in the vibration intensity we are aware of, 
was carried out by means of numerical simulations of a model for compaction || , and the data indicated the presence 
of memory effects in the evolution of the density of the system. Very recently M], memory effects were also directly 
' observed in a series of experiments. The results showed that the system has a short-term memory of its shaking 
. history, so that the response in the evolution of the density to a change in the vibration intensity at a given time, is 
not determined by the density at that time. Mathematically, this phenomenon implies that the time evolution of the 
density does not obey a closed ordinary first order differential equation. 

In this paper, we propose a general theoretical framework to understand the origin and characteristic features of 
the memory effects seen both in simulations and in experiments. Using quite plausible hypothesis, we will be able to 
explain the short-time response of the system to a small change in the vibration intensity. In particular, the theory 
predicts that a decrease (increase) in the intensity can lead to an increase (decrease) of the compaction rate on short 
time scales, in agreement with experiments. Nevertheless, this has not to be necessarily so. If the change in the 
intensity is made at the early stages of the compaction process, the theory we will develop leads to a modification of 
the compaction rate having the same sign as the intensity change. In fact, there is a critical time, which depends on 
the tapping intensity before the change, separating the regions of "normal" and "anomalous" responses. The existence 
of these two different regimes has not been verified experimentally up to now, although such a behaviour has been 
numerically observed in a simple model for granular compaction ||. 

As an illustration of the theory, we discuss its application to a model for compaction introduced recently . The 
model has already been showed to reproduce the qualitative behaviour of granular materials under tapping. Here we 
will show that it also captures the same short-term memory effects seen in the experiments. Moreover, it fits perfectly 
in the general scheme developed in this paper, therefore providing a first test of validation of the ideas in which the 
theory is based. We have also used the model to investigate the relaxation of the system following a perturbation in 
the vibration intensity for a short time period. This idea also originates from the experiments reported in ref. [0. 
The results indicate that the response function is accurately described by a Kohlrausch- Williams- Watts (KWW) or 
stretched exponential function. 

The paper is organized as follows. In the next Section, some general properties of the equation governing the time 
evolution of the density in tapping processes of granular media are discussed. These properties are used in Sec. [II to 
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analyze the short-term memory effects by considering the response of the s yst em to a small change in the vibration 



intensity. The theory is particularized for a simple model for tapping in Sec. IV, where other patterns of change of the 
vibration intensity are also considered. The choices were originated from the experiments reported in r ef. g). Finally, 



1 



Sec. contains some additional comments and final remarks, as well as a relevant discussion of the relationship of 
our work to previous experimental and theoretical studies. 

II. EVOLUTION OF THE DENSITY IN DISCRETE TAPPING PROCESSES 

Let r denote the dimensionless parameter characterizing the intensity of the vibration applied to the granular 
medium. In typical experiments Q], T is defined as the ratio of peak acceleration of a tap to the gravity g. Under 
very general conditions, the time evolution of the density p in a discrete tapping process will be given at a mesoscopic 
level by an equation of the form 

p^^ = / 1 (r) Ml (t)-/ 2 (r) / i 2 (t). (i) 

Here the time is measured in units of complete taps in some continuous limit, fi(T) and / 2 (r) are semidefinite positive 
functions of T having dimensions of frequency, and p\ (t) and /x 2 {t) are positive quantities depending on the state of 
the system, but they are not univocally determined by the density at the same instant p(t). Therefore, Eq. (|l|) is not 
in general a closed equation and can not be solved by itself. The two terms on the right hand side of the equation 
describe elementary processes increasing and decreasing the density, respectively. 

The structure of Eq. ([[]) as a gain-loss equation is consistent with the experimental observations in compaction 
processes, as we will discuss in detail in the following. Also, if the elementary events taking place in the system being 
vibrated can be described by means of a master equation, a formal equation like this follows directly. This is the case 
for some simple kinetic models for compaction introduced recently ■ 

Since Eq. (Q) describes the evolution of the density as a consequence of tapping, the functions /1 and / 2 must vanish 
in the limit of no tapping T = 0, so that 

/i(0)=/ 2 (0)=0. (2) 
Because of continuity, it follows, at least for small values of the intensity T, that 

/{(T) = A/ 1 (r)>0 l &(T) = A /a (r) > 0. (3) 

We will assume that the above inequalities hold for arbitrary T. The physical reason for this assumption is that we 
expect the number of elementary processes taking place in the system to increase as T increases. Of course, this does 
not imply by itself that the rate of variation of the density also increases. The behaviour of p depends on the net 
balance between the gain and loss elementary events, as indicated by Eq. (Q). This picture is in agreement with the 
qualitative role of temperature played by the shaking intensity in many different aspects P,f7|[io|-|l"3]] . 

In the long time limit of a tapping process with constant T, the experiments show that the system reaches a steady 
state with a density p s , which is a monotonic decreasing function of T, as displayed by the "reversible" branch in 
cycling experiments [|9|. Let us point out that the relaxation process is very slow, and for very small values of T 
the steady density is hard to reach within the experimental time scale. Therefore, the function /O s (r) verifies that 
dp s (T)/dT < 0, and it is bounded by the two formal limits 

Pmin = lim Ps(r), p m ax = Mm /O s (T). (4) 
T— >oc r — >0 

Particularization of Eq. ([!]) for a steady state yields 

A(I> ls = / 2 (I> 2s , (5) 

where p\ s and p2s denote the steady values of the quantities p\ and P2, respectively. As pointed out above, H\{t) and 
/Lt 2 (£) are not expected to be simply functions of p(t) in general. But, on the other hand, it seems sensible to assume 
that the steady state reached by a given system in a tapping experiment is fully determined by the intensity T or, 
equivalently, by p s . Therefore, we assume that p± s and /x 2s are functions of p s , and in the following we are going to 
investigate some qualitative properties of these functions. For p s — > p m in, p2s{ps) must vanish, since by definition at 
Ps = pmin there are no processes decreasing the density. Therefore, it is 

lim P2s(Ps) = 0, lim pi s (p s ) > 0. (6) 
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The second relationship expresses that starting from the density p m im any tapping process of arbitrary intensity T 
can only produce an increase of the density. What happens in the steady high density limit?. A similar argument to 
the one carried out above would lead to 

lim p ls (p s ) = 0, lim p 2s {Ps) > 0- (7) 

Nevertheless, some care is required when analyzing this limit. Simple models for discrete tapping lead to an absorbent 
steady state in the high density limit That means that the system will not be able to leave this state when 

submitted to tapping of arbitrary intensity. This is equivalent to say that P2s(ps) also vanishes for p s — > p ma x- As a 
consequence, and in order to include such a possibility in our formulation, instead of Eq. (|7|) we will assume the more 
general and precise condition 

rim ^4=0, (8) 
i. e., p2s >> pis when p s — > p max , and the density loss term is dominant in that limit. Let us note that Eqs. (((]) yield 

hm - — - — oo. (9) 

p,->Pmlm P2s{ps) 

The simplest behaviour that is consistent with Eqs. (^) and (^) is that the ratio p\ s j P2s be a monotonic decreasing 
function of the steady density p s going from infinity to zero. Since there is not any physical reason to expect a more 
complicated density dependence, we assume this is the case in our formalism. From the steady condition given by 
Eq. (§), it follows that 

Pis /a(r) _ , , , s 

= TTF\ =9( T >- ( 10 ) 

P2s fi[T) 

The function g(T) is a measure of the rate of the decompaction processes with respect to the rate of the compaction 
ones. Because of Eqs. (||) and (||), it is 

lim g(T) = oo, lim g(T) = 0. (11) 
Taking derivative with respect to the intensity T in Eq. ( |Io| ) it is obtained: 

dl dl dp s \ p 2s J 

where we have taken into account the monotonically decreasing density dependence of p\ s / P2s assumed above. The 
physical meaning of Eq. ([l2]) is evident: the rate of the decompaction processes grows faster with T than the rate of 
the compaction processes. Also this implication of our assumptions seems physically plausible. 

In summary, we can write the equation for the time evolution of the density in discrete tapping processes as 

( ^ = h (T)[p l {t)-g{T)p2{t)], (13) 

with fi and g being positive increasing functions of T, both of them vanishing in the limit r — > 0. The quantities 
pi (t) and p2 (t) are some moments of the complete distribution function of the system, and they contain the influence 
of correlations on the evolution of the density. As a consequence, ( filf ) is not a closed equation. 

Because g(T) vanishes for T — > 0, if a tapping experiment with low enough intensity T is carried out, the decom- 
paction term /i(r)g r (r)/x 2 (t) will be negligible in the first stages of the process, i. e., 

Mt) »1, (14) 



S(I>2(0 

and the evolution of the system will be approximately described by 

dp(t) 



dt 



A(I>i(t). (15) 
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At much later times, when pit") is close enough to the asymptotic steady value, the decompaction contribution in Eq. 
( |l3| ) plays a decisive role, leading to a steady density p s < p ma x, and it is 

MlW =0(1). (16) 



s(i> 2 (t) 



The observed behaviour that the system tends towards a steady state and, therefore, a regime where p\ 1 p 2 , and g 
verify a relationship of the form given in Eq. (|l6|), can be understood if pi(t) decreases in time while p 2 increases. 
Quite interestingly, this is consistent with a mean field approximation in which p\{t) is replaced by pi s [p(t)] and 
t^2(t) by p2s[p(t)]. Since p(t) increases monotonically in time, and p\ s / P2s is a monotonic decreasing function of the 
density, it follows that the left hand side of Eq. (|l^) will decay in time. 

Of course, as long as Eq. (|l^) is accurate, the larger V the faster the compaction of the system, in agreement with 
experiments. Over a larger time scale, the complete Eq. ([l5|), including the decompaction term, is needed in order to 
explain the dependence of the steady density p s on T, and also the existence of a a slow long time tail in the relaxation 
of the density, once p s — p{t) is very small. In this context, the presence of an "anomalous" density relaxation, following 
an inverse logarithm law, would be associated to some specific dynamical properties of the compaction term p± when 
the system is submitted to "nonlinear" tapping processes jl],0|. We use the term "nonlinear" here in the sense that, 
in the experiments, the initial value of the density is not very close to the steady density. 

Later on, we will show that an evolution equation like Eq. (^) applies in the case of a simple model recently intro- 
duced to describe discrete tapping . Another similar equation is obtained for the "parking" model (|||[l4|,|l5| , 
although this latter refers to continuous vibration processes, in which the system is not allowed to relax to a metastablc 
configuration between every two vibration cycles. In the parking model, the state of maximum density p max is not 
totally absorbent, but this possibility has been included in our theory, as discussed below Eq. (^). By identifying the 
intensity of tapping T with the ratio between the desorption and adsorption rates in the parking model, it is trivial 
to check by using the expressions in Ref. j^] that the quantities corresponding to p\ s and p 2s verify that their ratio is 
an increasing function of p s , and also the limiting behaviour given in Eqs. (g) and (^|). Moreover, the steady density 
is an increasing function of the quantity playing the role of the vibration intensity. In conclusion, the parking model 
belong to the general class of systems we have considered. 



III. RESPONSE TO SMALL VIBRATION INTENSITY JUMPS 

In this Section we will investigate whether Eq. (|T3|), which has been built under very general arguments and is 
expected to have a wide range of applicability, is able to predict the memory effects recently observed in vibration- 
induced compaction in granular materials Q. The fact that the equation is not closed for the density, implies that 
its time evolution in a given experiment with constant T is not determined by its the initial value. Starting from the 
same value pa, different time evolutions are possible depending on the way in which the system was prepared. Our 
aim is to analyze some particular relevant manifestations of this general statement. 

Consider that, starting from a given configuration, the system is tapped with an intensity T. At a certain time t Wl 
the intensity is instantaneously changed to T + AT. Quite peculiarly, it has been observed in the experiments that 
the change in the compaction rate has opposite sign that Ar on short time scales, though in the long time regime the 
relaxation is slower for smaller values of the intensity of vibration T. The same kind of effect has also been previously 
found numerically in some models for compaction |5|,jl6LEI , although it only shows up when the time interval t w is not 
too short. If r is changed at the beginning of the compaction process, the variation of the compaction rate has the 
same sign as Ar |5| JT7|. 

Application of Eq. (|13|) for the instant t~, just before the change in the intensity of the vibration, yields 

T w = pit') = /i(r) [pl w - g(T)p 2w ] , (17) 

where p^ w = pi(t~) and p 2w = p 2 (t~). When the intensity of vibration is changed into T + AT, the compaction rate 
becomes 

r' w ee pit+) = /i(r + Ar) [p+ w - gir + Ar>+J • (is) 

The continuity of the distribution function of the system implies that /i7 = pi and p^ = p\ for an instantaneous 
jump of r, although there is a discontinuity Ar w = r' w — r w in the compaction rate. For Ar small we can approximate 

^ r w i i /'(r) , 

-j^r = fi( T ) i^iw - g(T)p2w] - .fi(T)g'(r)p 2w = JJj^ r ™ - fi( r )9'{T)v>2w (19) 
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Therefore, if over the compaction curve corresponding to intensity T we define the function 

X(t) = j^r(t)-f 1 (T)g'(T)p 2 (t), (20) 

the sign of this function at the time t w when the intensity is changed will determine the relative behaviour of Ar w 
with respect to Ar, for infinitesimal changes of the latter. If X w = X(t w ) < 0, the anomalous response observed in the 
experiments will follow, while if X w > the compaction rate will change in the same direction as Ar. Let us analyze 
the sign of the function X(t). In the long time limit, formally t w — > oo, the system is known to reach the asymptotic 
steady density, so that r w — > and, consequently, 

Ao. = lim x(t) = -A(rV(i> 2s (r) < o, (21) 

t — >oc 

where it has been taken into account that both /i(r) and £i(r) are positive increasing functions of T and that 

M 2s (r)>o. 

On the other hand, if the initial density in the experiment is the minimum possible density of the system at rest 
Pmin , corresponding to the random loose packing configuration, it follows from the properties of /i2 that 

}™m = TTT,r(0)>0, (22) 
in 1 J 

Even though we have considered in our discussion that p(t = 0) = p m in in order to derive the above inequality, the 
same result will apply if the initial density is close enough to it, so that the first term on the right hand side of Eq. 
( pp| ) dominates the second one in the initial regime. 

Then, we conclude that for short times Ar w and Ar have the same sign, while for large times their signs are 
opposite. This renders compatible and explains what is seen in the experiments and also in numerical studies of 
simple models. From our analysis it follows that there is (at least) a time t Cl which depends on the value of T, such 
that the response of the system to a small variation of the intensity of tapping is qualitatively different for t < t c and 
t > t c . 

The study carried out in this Section has been restricted to small instantaneous changes in T, allowing the use 
of a linear analysis of Eq. ( |l3| ) . Whether the behaviour of the system remains the same when submitted to a finite 
change in the shaking intensity, it can not be inferred from our analysis. In this case, nonlinear effects can modify 
dramatically the response of the system. More will be said about this in the next section of the paper. 



IV. APPLICATION TO A SIMPLE MODEL FOR COMPACTION 



The general scenario developed in the previous sections will be particularized here for a one-dimensional lattice 
model for compaction ffljio]]. In the model, each site i can be either empty or occupied by a particle. A variable mi 
is defined, being rrii = 1 in the former case and = in the latter. A configuration of the system is fully specified 
by giving the values of all the variables m = {m^}. As usual, we will refer to the empty sites as being occupied by a 
hole. 

Let us describe the dynamics of the system when submitted to a discrete tapping process. Mechanical stability 
requires that all the holes be isolated, i.e. surrounded by two particles, at the end of every tap. The time evolution 
of the system is defined as a Markov process, and formulated by means of a master equation for the probability 
distribution of the system 0,0. The equation contains the transition rates IT(m|m') from state m' to state m. 
There are three kinds of possible transitions. Indicating only the variables associated to the sites involved in the 
transitions, the nonvanishing transition rates are: 

1. Elementary diffusive events conserving the number of particles, 

TF(010|100) = VF(010|001) = |, (23) 

2. Transitions increasing the number of particles, 

W(010|101) = -, (24a) 

W(001|101) = PF(100|101) = (24b) 
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3. Transitions increasing the number of holes, i.e. decreasing the number of particles, 



VF(01010|00100) = —, (25a) 



W(01010|01000) = W(01010|00010) = — . 



(25b) 



In the above equations, a is a positive constant, characterizing the tapping process completely, and playing in the 
model a role similar to the intensity of vibration Y in real experiments. For a ^ 0, the system evolves from any 
arbitrary initial configuration to a final steady state with density 



p.(a) = i|l + (l + 4a) 1/2 



From here it follows that 



being 



lim p s = 1 = pmax, lim Ps = ~ = Pmin, 
a— >0 a^oo 2 



da 



(26) 



(27) 



(28) 



for all a. Therefore the density in the model has the same kind of dependence on the intensity a as assumed in the 
general discussion in Sec. The time evolution of p is obtained from the master equation for the model, and reads 



ax wl (t) 



a 

T 



2^00100 (t) + —xoxooo(t) 



1 



^oooio(^) 



(29) 



where xoio is the concentration of three-site clusters of the form hole-particle-hole, xooioo ^ s ^ ne concentration of 
five-site clusters formed by a hole between two pairs of particles, and so on. Comparison of Eqs. (|l|) and (|9|) allows 
to identify 



fi(a)=a, f 2 (a) 



(30) 



Pi{t) = xxox(t), p 2 (t) = ^ooioo (*) + 2 [^QioooW + »oooio(*)] ■ (31) 

In the steady state, the only correlations in the system are those forbidding to have two nearest-neighbour holes 
L8[. Then, it is a simple matter to compute the steady values of the several cluster concentrations appearing in Eq. 
3l| ) with the result 

= (32) 



In the limit p s — > p m in = 1/2, 



pi 



Mis -> 2' ' i2s ~* °' ( 34 ) 



while in the high density limit p s — > p m ax = 1, both p\ s and p 2s vanish, as a consequence of the absorbent character 
of the state with all the sites occupied by particles. The ratio 



G 



m M (2 Ps -iy { ) 

vanishes in this latter limit. Equations ([mJ) and ( |35| ) are in agreement with Eqs. (0) and (||). Moreover, p-is(ps) / P-2s(Ps) 
is a monotonic decreasing function of p s and, consistently (see Eq. (p^Oj)), 

is an increasing function of a, vanishing in the limit a — > 0. 

We conclude that this model for compaction fits perfectly the general picture developed in the previous Sections. 
Equation ( |l3|) particularized for the model is 

—- = a(m - afja), (37) 

with /Lti and \xi defined in Eq. (|3l|). To solve Eq. ( |37| ) we would need some (approximate) expressions for the cluster 
concentrations as functions of the density. 

If we submit the system to the tapping experiment described in Sec. |lTl| , the effect of the intensity change Act at 
t = t w on the compact rate will be given by 

= — -ap 2w , (38) 



Ac 



in the limit of small Aa. Therefore, the function determining whether the response of the system will be "normal" 
or "anomalous" is 

X(t) = ^--a t x 2 (t). (39) 
a 

In Fig. [l] this function is plotted for a = 0.15. The curve has been obtained by Monte Carlo simulation of the 
master equation of the system. The data represent an average over 10 different runs. The initial state was the one 
corresponding to the steady minimum density. For this particular value of the intensity a, X(t) changes sign between 
taps 19 and 20, i.e. 19 < t c < 20. For comparison purposes, we have also plotted the mean field approximation 
for A(i) (dashed line). The latter has been constructed by substituting in Eq. ( |39| ) and liiit) by /ii s (p(t)) and 

p2s(p(t)), respectively, and using for the density the simulation results. It is seen that the mean field approximation 
also changes sign, but for larger times, and it is always above the "exact" Monte Carlo curve. This is consistent with 
the mean field approximation giving a faster approach to the steady state than the actual relaxation of the system 

& 

According to the results derived in this paper, Ar w is expected to have different sign that Aq for t w > t c and the 
same for t w < t c . In order to check this theoretical prediction, we have carried out series of Monte Carlo simulations, 
all them starting in the minimum density configuration, with a — 0.15. At t w = 50, the value of the intensity 
a was instantaneously changed to a'. The results for four different values of a' are reported in Fig. |^, namely 
0.1,0.125,0.15,0.175, and 0.2, from top to bottom. The central value corresponds to no change. Since in these 
simulations it is t w > t c , the compaction rate is observed to decrease as the value of a' increases. It is also seen that 
the amplitude of the jump in the compaction rate is larger for a' = 0.2 (Aq = 0.05) than for a' = 0.1 (Aa = —0.05). 



This feature can not be explained by Eq. (38), and it is due to nonlinear effects that have been neglected in the linear 
approximation used here. This will be analyzed below. 

In Fig. H the same set of experiments is carried out, with the only difference that in this case the intensity a is 
modified at t w = 10 < t c . The several curves correspond to the same values as in Fig. pi but now they are ordered 
from bottom to top. As predicted by the theory, the variation of the compaction rate has the same sign as the change 
in a. Moreover, the same kind of nonlinear effects as in Fig. ^| are present. 

Now, we will briefly discuss the nonlinear corrections in Aa to the change in the compaction rate. It is easy to 
show that 

Ar w = Aa\ w - (Aa) 2 p 2w . (40) 



The second term on the rhs of Eq. (40) is neglected in the linear approximation. In this simple model, the nonlinear 
correction is always negative, so it can modify dramatically the response of the system to the jump Aa if the linear 
term AaA^ > 0. In particular, there is a critical value 
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Plw 

such that Ar w = 0. For smaller jumps, |Aa| < |Aa c |, the sign of Ar w is the one predicted by the linear approximation, 
but for larger jumps, |Aa| > |Aa c |, the sign of Ar w is the opposite to the prediction of the linear approximation. 
For the sake of concreteness, in Fig. |] we have repeated the numerical experiment of Fig. ||, but with larger intensity 
jumps. From the Monte Carlo simulation, we obtain the critical value Aa c ~ —0.1 for a = 0.15 and t w = 50. Then, 
at — — 50 we change the vibration intensity from a = 0.15 to a' = a + Aa c = —0.05, finding that the compaction 
rate does not change in the short time limit t — t w — ► 0. Moreover, if the vibration intensity is further decreased, 
a' = 0.03, the compaction rate also decreases, while the linear approximation predicted an increase of the compaction 
rate if a' < a, since X w < 0. 

Following Ref. Q, we have also considered another series of numerical experiments where the system was tapped 
up to the same density with three different intensities, a — 0.2, 0.15, and 0.1, respectively. Afterwards, the system 
was always tapped with the same intensity a' — 0.15. The time evolution of the density is shown in Fig. [5], where the 
time origin for each experiment has been taken at the time when the system reached the prescribed density, namely 
p = 0.8. The figure clearly shows that the evolution of the density for t > strongly depends on the previous tapping 
history, indicating the relevance of short-term memory effects. Mathematically, this is equivalent to say that pi(t) and 
(J,2(t) in Eq. ([l]) are not determined univocally by the density at the same time, so that it is not in fact a closed first 
order ordinary differential equation. Note that in all the plotted curves the jump in the compaction rate has opposite 
sign than the variation of the intensity. We have verified that X(t) is negative at the time in which the intensity is 
modified in all cases, the behaviour being then consistent with the theory. 



V. DISCUSSION 



Along this paper, we have studied the non-equilibrium linear response of a vibrated granular system to an instan- 
taneous change in the intensity of the taps. In the first part, a general theory was developed on the basis of some 
plausible hypothesis about the mesoscopic dynamics of the system. The results are in qualitative agreement with the 
experimental observations. In particular, the presence of short-term memory effects appears as correlated with the 
relaxation properties of the system at constant intensity. An important theoretical prediction, not observed in the 
experiments yet, is the existence of a critical time t c . For times t < t c the response of the system to a change in the 
intensity is "normal" , in the sense that an increase in the intensity produces a positive jump in the compaction rate. 
On the other hand, for t > t c an "anomalous" response is produced. The change in the compaction rate has opposite 
sign that the modification of the vibration intensity, in contrast with the long time behaviour found in experiments, 
where the relaxation is faster for larger vibration intensity. 

In the second part of the paper, a simple model for compaction has been considered. It is shown to fit perfectly into 
the general scheme developed before, allowing a detailed quantitative analysis of the theoretical predictions. This is 
not a peculiarity of this model, since the "parking" model || also verifies all the conditions assumed in the theoretical 
framework. In fact, this is not surprising because this latter model has a mathematical structure very similar to the 
one considered in this paper ||. 



In Sec. IV we have shown that our model reproduces the experimentally observed behaviour of the system when 
submitted to changes in the vibration intensity under different conditions [|4[ . Now we will refer to a more complicated 
pattern of changes in the intensity that are also discussed in ref. Q|. First, the system is shaken with an intensity 
Tq (ao in the model notation) for a long period of time, so that the system practically reaches a steady density. 
Afterwards, at a time taken as t = 0, the intensity is switched to T > r for a given period of time t and, finally, 
the system is tapped again with the original intensity To, and the subsequent relaxation of the system is studied. 
Experimentally it was found that the relaxation is slower the larger to, the system "ages". Moreover, on the basis of 
a simple two-state model, it was proposed that 

exp(-Ko^) 

P{t)-Ps , (42) 

for t — to 3> 1. In the above expression kq is a decreasing function of to- Josserand et al. Q also reported that the 
relaxation can be fitted by a superposition of exponentials, all of them with the same amplitude. We have carried out 
numerically this kind of experiments in our model. In Fig. ^ we present the results obtained with a = 0.3, a — 0.5, 
and four different values of to, namely to — 1,2,4, and 8 from bottom to top. The plotted response function <p(t) is 
defined as 

Ps{a ) - p(t ) 
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where p(to) is the density of the system at the time in which the the intensity is switched back to ao- For the model, 
the steady values of the density can be computed analytically ||, and it is p s {a — 0.3) ~ 0.8371. For a given time 
interval t — to, <p(t) increases with the "waiting" time to. Thus, the relaxation is slower for larger to, consistently 
with the experimental observation M. Also, we have fitted (solid lines) the data to a stretched exponential or KWW 
function [HI, 



PKWw{t) = p s (ao) - [Ps(«o) - p(to)] exp 



t-t 



(44) 



with r and [3 being fitting parameters. As observed in the figure, the fit is quite satisfactory, except for times very 
close to to and, probably, for very large times. The parameter /3 in Eq. (Q) measures the width of the relaxation time 
distribution, he values we have found go from (3 — 0.366 for to = 1 to (3 = 0.478 for to — 8. The latter is close to the 
value 1/2, characteristic of systems whose dynamics is dominated by one-dimensional diffusive processes. However, 
the KWW relaxation is not equivalent to a superposition of exponentials with the same amplitude, as proposed in [Q. 
Therefore, this point deserves more work in the future, both theoretically and experimentally. With respect to the 
long time behaviour predicted by Eq. (^) we could not reach a definite answer. Although the numerical data seem 
to be compatible with it, the noise is too large and further high-precision studies would be required. 

Finally, a crucial point in the analysis presented in this paper is the small amplitude of the perturbation in the 
vibration intensity Ar. As pointed out at the end of Sec. HI, the behaviour following a large change in the intensity 
may be different. In the model considered in Sec. IV the nonlinear corrections are very simple, leading always to a 



decrease in the compaction rate and to the appearance of a critical value of the intensity jump, such that no change 
in the vibration intensity is observed in the short time regime. Moreover, for jumps larger than the critical one, the 
sign of the change in the compaction rate is reversed as compared with the prediction of the linear approximation. 
We think that it is worth looking for this kind of behaviour in other models for compaction, and also in experiments 
with real granular systems. 
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FIG. 1. Time evolution of the function A defined in Eq. ^jjj, for a vibration intensity a — 0.15 (solid line). Also plotted is 
the mean field approximation for A (dashed line). 

FIG. 2. Evolution of the density, p, as a function of the number of taps, t. Five numerical experiments are shown, the 
vibration intensity a was changed at t w = 50, where A < 0, from 0.15 to 0.1, 0.125, 0.15, 0.175, and 0.2, from top to bottom. 
Thus, the central curve corresponds to no change in the tapping intensity a (solid line), while the dotted and dashed lines 
correspond to a decrease and an increase in a, respectively. The "anomalous" response experimentally observed shows up. 

FIG. 3. The same experiment of Fig. |^, but the change in intensity is introduced at an earlier time t w = 10, at which 
A > 0. The curves correspond to the same values of a' as in Fig. but now are ordered from bottom to top. In this region 
the response is "normal", i. e., the compaction rate increases with the vibration intensity. 

FIG. 4. Time evolution of the density p when the vibration intensity a is changed at t w = 50, where A < 0, from q=0.15 
to a' — 0.05 (circles) and 0.03 (squares), respectively. The curve corresponding to a constant vibration intensity a = 0.15 is 
plotted for reference (solid line). For such large jumps, the linear approximation is not valid, and the compaction rate do not 
increase. In fact, a' — 0.05 corresponds to the critical value Aa c , for which no change in the compaction rate is observed for 
short times. 

FIG. 5. Time evolution of the density for a system which was tapped up to the same density p = 0.8 using three different 
intensities, a = 0.2 (circles), 0.15 (squares), and 0.1 (triangles). Afterwards, the system was always tapped with a = 0.15. The 
time origin for each experiment has been taken at the time when the system reached the prescribed density, namely p = 0.8. 
The evolution for t > strongly depends on the pre-history of the system. 

FIG. 6. Relaxation function (f>(t), defined in Eq. (^), of the model, when it is prepared by tapping for a long time with 
a = 0.3, and afterwards tapped for to = 1, 2, 4, and 8 (from bottom to top) with a larger intensity a' = 0.5. Finally, the 
intensity is turned back to the original intensity a — 0.3. All the curves tend to zero in the infinite time limit, and the solid 
lines are the best numerical fits to a KWW function. 
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